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Abstract 


In the present paper we propose a new MCMC algorithm for sampling from the posterior 
distribution of hidden trajectory of a Markov jump process. Our algorithm is based on the 
idea of exploiting virtual jumps, introduced by Rao and Teh (2013). The main novelty is 
that our algorithm uses particle Gibbs with ancestor sampling (PGAS, see Andrieu et al. 
(2010); Lindsten et al. (2014b)) to update the skeleton, while Rao and Teh use forward 
filtering backward sampling (FFBS). In contrast to previous methods our algorithm can 
be implemented even if the state space is infinite. In addition, the cost of a single step of 
the proposed algorithm does not depend on the size of the state space. The computational 
cost of our methood is of order 0(NE(n)), where N is the number of particles used in the 
PGAS algorithm and E(n) is the expected number of jumps (together with virtual ones). 

The cost of the algorithm of Rao and Teh is of order 0(\ Aj 2 E(n)), where \X\ is the size of 
the state space. Simulation results show that our algorithm with PGAS converges slightly 
slower than the algorithm with FFBS, if the size of the state space is not big. However, if 
the size of the state space increases, the proposed method outperforms existing ones. We 
give special attention to a hierarchical version of our algorithm which can be applied to 
continuous time Bayesian networks (CTBNs). 

Keywords: Continuous time Markov processes, Bayesian networks, MCMC, Sequential 
Monte Carlo, Hidden Markov models, Posterior sampling, CTBN 

1. Introduction 

Markov jump processes (MJP) are natural extension of Markov chains to continuous time. 
They are widely applied in modelling of the phenomena of chemical, biological, economic 
and other sciences. An important class of MJP are continuous time Bayesian networks 
(CTBN) introduced by Schweder (1970) under the name of composable Markov chains and 
then reinvented by Nodelman et al. (2002a) under the current name. Roughly, a CTBN is a 
multivariate MJP in which the dependence structure between coordinates can be described 
by a graph. Such a graphical representation allows for decomposing a large intensity matrix 
into smaller conditional intensity matrices. 

In many applications it is necessary to consider a situation where the trajectory of a 
Markov jump process is not observed directly, only partial and noisy observations are avail¬ 
able. Typically, the posterior distribution over trajectories is then analytically intractable. 
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The present paper is devoted to MCMC methods for sampling from the posterior in such a 
situation. 

In the literature there exist several approaches to the above mentioned problem: based 
on sampling (Boys et ah, 2008; El-Hay et al., 2008; Fan and Shelton, 2008; Fearnhead and 
Sherlock, 2006; Hobolth and Stone, 2009; Nodelman et ah, 2002b; Rao and Teh, 2013, 2012; 
Miasojedow et ah, 2014), based on numerical approximations (Cohn et ah, 2010; Nodelman 
et ah, 2002a, 2005; Opper and Sanguinetti, 2008). Some of these methods are inefficient, like 
modification of likelihood weighting (Nodelman et ah, 2002b). Other approaches involve 
expensive computations like matrix exponentiation, spectral decomposition of matrices, 
finding roots of equations. There are also approximate algorithms based on time discretiza¬ 
tion. To the best of our knowledge the most general, efficient and exact method is that 
proposed by Rao and Teh (2013), and extended to a more general class of continuous time 
discrete systems in Rao and Teh (2012). Their algorithm is based on introducing so-called 
virtual jumps and a thinning procedure for Poisson processes. In our approach we combine 
this method with particle MCMC discovered by Andrieu et ah (2010). More precisely, in¬ 
stead of forward filtering backward sampling algorithm used in the original version, we use 
particle Gibbs (Andrieu et ah, 2010) with added ancestor resampling proposed in Lindsten 
et ah (2012, 2014b). The proposed method is computationally less expensive. Moreover, 
our algorithm can be directly applied when the state space is infinite, in opposition to Rao 
and Teh (2012, 2013). 

2. Markov jump processes 

Consider a continuous time stochastic process {X(t),t> 0} defined on a probability space 
(17, P, P) with a discrete state space X. Assume the process is time-homogeneous Markov 
with transition probabilities 

P\s, s') = P (X(t + u) = s'\X{u) = s) , 

for s, s' E X. The initial distribution is denoted by v(s) = P(X(0) = s). Since X is discrete, 
v can be viewed as a vector and P t as a matrix (both possibly infinite). The intensity matrix 
is defined as follows 

Q(s,s') = lim^[P*(s,s') -1(s, s')] , 

where I = P° is the identity matrix. Equivalently, Q(s, s') is the intensity of jumps from s 
to s', i.e. 


P (X(t + d t) = s / |A(t) = s) = Q(s, s')dt for s ^ s' , 

P (X(t + d t) = s\X(t) = s) = 1 - Q(s)dt , 

where Q(s ) = — Q(s,s ) = ^2 S /^ S Q{s, s') denotes the intensity of leaving state s. Clearly 
we have ^2 S , Q(s,s') = 0 for all s E X. Throughout this paper we assume that Q is non¬ 
explosive (Norris, 1998), which means that almost surely only finite number of jumps occur 
in any bounded time interval [0,f max ]. This assumption is fulfilled in most applications we 
have in mind. Sufficient and necessary conditions for Q to be non-explosive can be found in 
(Norris, 1998)[Thm. 2.7.1, 2.7.2, Cor. 2.7.3]. From now on the interval [0,t max ] is fixed. Let 
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there be m jumps and let these jumps occur at ordered moments T = (t\,, t m ). Moments 
of jumps T with a corresponding sequence of states, denoted by S = (so, si, • .., s m ) = 
(X(0), X(t\),..., X(t m )), fully describe the trajectory X([0, t max ])- By definition of the 
process {X (t)}, every interval between jumps, tj—tj- 1 for j = 1 ,..., m, with the convention 
to = 0, is distributed according to the exponential distribution with parameter Q(sj- 1). The 
skeleton S' is a discrete time Markov chain with initial distribution v and transition matrix 
given by 


Q(s, -s') 
Q(s) 


if s 7 ^ s' ; 
if s = s' . 


Thus random variable (T, S) has density 


D( A 

p(T,S) = u(s 0 ) Q(sj-i)exp{-Q(sj)(tj - tj- 1 )} eX p{-(5(s m )(t max - t m )} 

j = 1 

m 

— r'(so) 11 Q(sj —l j Sj ) exp { Q(sj ) (tj tj—i)} exp { Q(s m )(t max tm)} j (1) 

3=1 


where rn = |T|. The last factor exp {— Q(s m )(t mSLX — t m )} comes from the fact that the 
waiting time for jump m + 1 can be arbitrary but greater than t max — t m . 


3. Virtual jumps 

Let {V(t)| be a homogeneous Markov process with intensity matrix Q and let R(s ) > Q(s ) 
for all s 6 X. Consider the following sampling scheme (Rao and Teh, 2012), based on 
dependent thinning, i.e. rejection sampling for an inhomogeneous Poisson process (Lewis 
and Shedler, 1979). We generate a sequence of potential times of jumps (t\,t 2 , ■ ■ ■)■ For a 
given moment tk-i and a current value of the process §k-i = X(tk~ i), we draw the next 
time interval t^ — tk-i from the exponential distribution with parameter R(sk-i). With 
probability Q(sk-i)/R(sk-i) the process jumps at time tf. to another state, and this new 
state is §k with probability Q(sk-i, Sk)/Q(sk-i). With probability (1 — Q(sk-i))/R(sk-i) 
the process does not jump and we put Sk = Sfc-i- The resulting redundant skeleton S = 
(§ 0 ) h, S 2 , • • ■) is therefore a Markov chain with transition matrix P defined by 


P(s,s') 


' Q(s,s') 
R(s) 

T Q(s) 

R(s ) 


if s ^ s' ; 


if s = s' . 


( 2 ) 


We summarize this procedure as the following algorithm 1. 
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Algorithm 1 Thinning procedure. 
Set to = 0 and k = 0. 

Draw so ~ u {')- 
while t k < imax do 

Set k = k + 1. 

Draw W ~ Exp(R(sk-\)) ■ 

Set tk = tk- 1 + W. 

Draw s k ~ P(h-i, • )• 

end while 


As before, consider the process in a fixed interval of time [0, U,]. Let T = (ti, £ 2 , ■ • ■, t n ) 

be the set of moments generated by the above algorithm. Let J = {k > 1 : Sk ^ Sfc-i}. 
Denote by T = Tj and V = T-j the moments of true jumps and virtual jumps, respectively. 
The process {X(t)} resulting from the algorithm has the same probability distribution as 
that in Section 2. This fact is explicitly formulated in Proposition 1 below. 

Proposition 1 The marginal distribution of (T = Tj,Sj ) has the density given by (1). 

This proposition is known and can be found e.g. in Rao and Teh (2012). However, to make 
the paper self-contained we give a proof in the Appendix. The following corollary shows 
that, conditionally on (T, Sj ), i.e. on true jumps and the skeleton, the set of virtual jumps 

V is an inhomogeneous Poisson process with intensity R(X(t)) — Q((X(t)). 

Corollary 2 Let Vj denote moments of virtual jumps between two adjacent true jumps tj -1 
and tj. The conditional density of Vj is given by 

p(Vj\X(tj-i) = s , tj—\ , tj') = (. R(s ) - Q(s))\ v i\ exp {-(tj - tj-i)(R(s) - Q(s))} . 

The proof is also given in the Appendix. In the sequel we will work with the redundant 
representation of the process {X(t)} introduced in this section. For simplicity, let us slightly 
abuse notation and from now on write S instead of S. Clearly, the density of (T, V, S ) is 
given by 


P(T, V, S ) = i/(s 0 ) R(sk-i)P(s k -i,s k ) exp 

k=1 
n 

= u(s 0 ) Il(fl(sn) - Q(sfc-i)) 1(sfc=Sfe - l) * * * V * * * * X g(s fc _i,Sfe) 1(sfc ^ fc - l) (3) 

k=1 
n 

X n exp {-R(s k )(t k - tk- 1 )} exp {-R(s n )(imax - t n )} , 
k= i 

where n = \T\ + |V| is the total number of jumps. 

Now we describe two particular choices of intensity R. The first one leads to the so- 
called uniformization (Jensen, 1953; Cinlar, 1975; Hobolth and Stone, 2009; Rao and Teh, 
2013). Let A > max s Q(s) and consider moments of potential jumps distributed according to 
homogeneous Poisson process with intensity A. Precisely, we consider the thinning procedure 


R(X(u))du 
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with R(s) = A. Clearly, the joint distribution of true jumps, virtual jumps and skeleton is 
now given by 


p(T, V, S) oc A'V(s 0 ) P(s k - 1 , s k ) , 


(4) 


k=1 


where P is the transition matrix of discrete time Markov chain defined by 

( Q(s,s') 

P{s,s') = l 


A 


if s / s' ; 


(5) 


1 - 


Q(s) 

A 


if s = s' 


In the case of uniformization, conditionally on the trajectory Al ([0, imax]); the virtual jumps 
form a piecewise homogeneous Poisson process with the intensity constant and equal to 
A — Q(X(tj )) on every time interval [tj-tj+i) for j = 0,1,..., |Tj. 

The second natural choice is to make virtual jumps distributed as homogeneous Poisson 
process with intensity 9 > 0. Let R(s) = Q{s) + 9. Then the thinning procedure leads to 
the following probability distribution: 


p(T, V,S) oc i'(so) P(s k -i,s k )(9 + Q(sk-i)) exp {-(Q(s fc _i) + 6)(t k - t k - 1 )} 

fc=i 

x exp {-(Q(s„) + 9){t max in )} (6) 

n 

= v(so) n Q( s k- 1 > Sfc) 1(sfc ^ sfc-l) ex P {-Q(s k -i)(t k - 4- 1 )} 

k=1 

x exp {-Q(s n )(t max in )}0 |y I exp {-9t max} 5 


where the transition matrix P of discrete time Markov chain which generates the skeleton 
S is given by 


' Q(s,s') 
Q(s ) + 9 


if s s' ; 


pM = < 


(7) 


9 

. Q( s ) + 9 


if s = s' . 


4. Continuous time Bayesian networks 

Let (V,£) denote a directed graph with possible cycles. We write w u instead of (w, u ) £ 
£. For every node w £ V consider a corresponding space X w of possible states. Assume that 
each space X w discrete. We consider a continuous time stochastic process on the product 
space X = EUv Xw- Thus a state s £ X is a configuration s = (s w ) = (s w ) w< =\;, where 
s w £ X w . If W C V then we write syy = (s m ) m gw for configuration s restricted to nodes in 
W. We also use notation Xyy = n t „ew so we can wr if e s w £ Xyy. The set V\{iu} 
will be denoted simply by — w. We define the set of parents of node w by 

pa(u;) = {u £ V : u -* tc} , 
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and we define the set of children of node w by 


ch(iy) = {u G V 

Suppose we have a family of functions Q h 


w —y w} . 

x Mw) X (Xw X X w ) -A [0, oo). For fixed 


c G <T pa ( w ), we consider Q w (c; •, • ) as a conditional intensity matrix (CIM) at node w (only 
off-diagonal elements of this matrix have to be specified, the diagonal ones are irrelevant). 
The state of a CTBN at time t is a random element X(t) of the space X of configurations. 
Let X w (t) denote its rcth coordinate. The process {X w (t) w£ \;,t > 0} is assumed to be 
Markov and its evolution can be described informally as follows. Transitions at node w 
depend on the current configuration of the parent nodes. If the state of some parent 
changes, then node w switches to other transition probabilities. If s w / s' w then 

If* ( X w (t T dt) — S w | X— w ( t) — S— w , X w {t ) — S w ) — Qw (Sp a (u!)) &w j S’w') 

Formally, CTBN is a MPJ with transition intensities given by 


Q(s,s') = 


Qw(s pa ( w ),s w , s' w ) if s- w = s'_ w and s w / s' w for some w, 
0 if S- w / s'_ w for all w, 


for s / s' (of course, Q(s, s ) must be defined “by subtraction” to ensure Q(s, s') = 0). 

For a CTBN, the density of sample path X = X([0, t max ]) in a bounded time interval 
[0, fmax] decomposes as follows: 


p(X) = v(X(0))l[p(X w \\X Mw) ) , 

w£V 


( 8 ) 


where v is the initial distribution on X and p(X w ||X pa (^) is the density of piecewise ho¬ 
mogeneous Markov jump process with intensity matrix equal to Q w (c, - , • ) in every time 
sub-interval such that X pa (^) = c. Formulas for the density of CTBN appear e.g. in (Nodel- 
man et al., 2002b, Sec. 3.1), (Fan et al., 2010, Eq. 2), (Fan and Shelton, 2008, Eq. 1) and 
(Miasojedow et ah, 2014). These formulas give a factorization of the main part of the 
density as in (8), but in the first three of the cited papers the initial distribution v is disre¬ 
garded. There are some subtle problems related to u, discussed in Miasojedow et ah (2014). 
Our notation p(X w \\X p ^ w ^) is consistent with the notion of “conditioning by intervention”, 
see e.g. (Lauritzen, 2001). Indeed, p(X w \\X pa («,)) is the density of the process X w at note 
w under the assumption that the sample paths at the parent nodes X pa f w \ are fixed and 
Xu,(0) is given, see e.g. Miasojedow et ah (2014), for details. Below we explicitly write an 
expression for p(X w \\X pa r w \) in terms of moments of jumps and the skeleton of the pro¬ 
cess (X Wl X Mw) ), as in (1). Let T w = (tg>... ,t ?,...) and T pa ^ = (t^,... ,*f (u °,...) 
denote moments of jumps at node w G V and at parent nodes, respectively. By conven¬ 


tion put t™ = t 


= 0 and t^ TW | +1 = = tmax- Analogously, S w and S pa ^ 

denote the corresponding skeletons. Thus we divide the time interval [0, t max ] into segments 
. _ 0; i ) ... |TP a ( w )| such that X pau ^ is constant and X w is homogeneous 
in each segment. Next we define sets Ij = {i > 0 : t pa ^’^ < tf < } with notation 
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JbegjJend for the first and the last element of Ij. Analogously to (1), we obtain the following 
formula: 


p(x w ||Xp,(„) )=p(T w ,S™ II S p,( ”>, ) 


n 


j =o 


n 

Aeij 


s 


w 

i—li 


sT) 


x JJ exp |-(if 

'beg} 


(9) 


x exp 




W +*■ 

Jbeg J 




pa(iu) # 


Formula (9) is equivalent to (Nodelman et ah, 2002b, Eq. 2) and (Fan and Shelton, 2008, 
Eq. 1), but expressed in terms of ( S,T ). 


5. Hidden Markov models 

Let X = {X(f), 0 <t< f max } be a Markov jump process. Suppose that process X cannot be 
directly observed but we can observe some random quantity Y with probability distribution 
L(Y|A([0,f max ])- Let us say Y is the evidence and L is the likelihood. We assume that the 
likelihood depends on X only through the actual sample path X([0, f max ]) (does not depend 
on virtual jumps). The problem is to restore the hidden trajectory of X given Y. From the 
Bayesian perspective, the goal is to compute/approximate the posterior 

p(X([ 0, t max ])\Y) oc p(X([0, t max ]))L(Y\X([0, f max ]))- 

Function L, transition probabilities Q and initial distribution v are assumed to be known. 
To get explicit form of posterior distribution we will consider two typical forms of noisy ob¬ 
servation. In the first model, the trajectory of X is observed independently at deterministic 
time points t* with corresponding likelihood functions L\(y\\X (t*)),..., Li(yi\X(t *)). 

Since X(t) is constant between jumps, for all* = 1,..., l we have Lj(yj| X{t*)) = L;(?/i|sj*), 
where i* = supjj : tj < t*}. If the trajectory of X is represented by moments of true jumps 
T, virtual jumps V and skeleton S then the posterior distribution is given by 

i 

p(T, V, S|Y) oc p(T, V, S ) H LiiviM (10) 

i= 1 

where p(T, V, S ) is given by (3). In Section 3 we considered two scenarios of adding virtual 
jumps: via uniformization and via a homogeneous Poisson process. The corresponding 
densities are given by (4) and (6), respectively. Observe that for both variants of adding 
virtual jumps, the posterior can be expressed in the following form: 

n 

p{S\T, V, Y) = is(s 0 )go(so ) P{s k -i, s k )gk(s k ), (11) 

k =1 

where P is a Markov transition matrix, g k are some functions which can depend on T, V, Y 
and n = |T| + \ V\. Hence the skeleton, conditionally on all the jumps and the evidence, can 
be treated as a hidden Markov model with discrete time. 
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In the second typical model of observation, the evidence Y is a fully observed continuous 
time stochastic process depending on X. Expressly, we assume that Y, given the trajectory 
of X, is a piecewise homogeneous Markov jump process such that the pair ( X,Y) is a 
CTBN with the graph structure X — > Y. Thus the likelihood can be expressed by (9), with 
X w = Y and X pa ^ = X. In this case we can easily obtain the same conclusion as before: 
for both variants of adding virtual jumps, the posterior can be expressed in the form ( 11 ). 
The skeleton is conditionally a hidden Markov model. 

6. MCMC algorithm 

Let us recall the standing assumption that evidence Y depends only on the trajectory of X 
but not on virtual jumps. This assumption covers most of usual scenarios and clearly implies 
that p(V\T, S,Y) = p(V\T,S). Now we are able to state the main algorithm. Similarly to 
(Rao and Teh, 2012, 2013), a single step of the iterative procedure is the following. We 
take the trajectory obtained in the previous step, represented by (T, S) (ignoring the virtual 
jumps). First we sample a new set of virtual jumps V. We can use two variants of sampling. 
In the case of uniformization, V is a piecewise homogeneous Poisson process with intensities 
A — Q(X(tk))- Alternatively, V is a homogeneous Poisson process with intensity 9. Next we 
generate a new skeleton S' using Markov kernel with p(S\T, V,Y) as invariant distribution. 
In (Rao and Teh, 2012, 2013), the authors use independent sampling by the forward filtering 
- backward sampling (FFBS) algorithm. In our approach we propose to use the particle 
Gibbs algorithm invented by Andrieu et al. (2010), which is described in the next section. 
Note that a new skeleton with fixed times of potential jumps leads to a new allocation of 
true and virtual jumps. So we obtain a new trajectory described by (T',V',S') such that 
Tub = T' U V'. Finally we remove virtual jumps to obtain a new state (' T',S'). The 
algorithm is summarized below. 


Algorithm 2 Single step of MCMC algorithm. 

Input: Previous state (T, S) and observation Y. 

1. Add virtual jumps V. 

2. Draw new skeleton S' from Markov kernel K(S, ■ ) targeting p(S\T, V, Y). Skeleton S' 
defines new allocation of virtual and true jumps T ', V' such T UV = T' U V'. 

3. Remove virtual jumps V'. 
return new state (T 1 , S'). 


The next proposition shows that this algorithm is ergodic. 

Proposition 3 Assume that A > max s Q{s) in the case of uniformization or 6 > 0 in 
the case of homogeneous virtual jumps. Assume that kernel K(S,S') leaves distribution 
p(S\T,V,Y) invariant and K(S,S') > 0 for all S' such that p(S'\T,V,Y) > 0. Then 
MCMC algorithm described above is f-irreducible, aperiodic with stationary distribution 
7 r(T, S) = p(T , S\Y). Thus for n-almost all initial positions, the algorithm is ergodic, i.e. 

\\M((T, S), ■ ) m — 7r(-)||t v ”•—>° 0 , 

where M denotes the kernel of our MCMC algorithm. 





Particle Gibbs for MJPs 


Proof By construction it is clear that p(T,S\Y) is a stationary distribution. Since 
K(S,S) > 0, with positive probability it happens that the skeleton does not change and 
hence virtual and true jumps remain unchanged. Clearly M((T, S), (T, S)) > 0 and so 
Markov chain M is aperiodic. The assumption A > max s Q(s) or 9 > 0 ensures that the 
step of adding virtual jumps can reach any configuration of virtual jumps. Together with 
the assumption K(S, S') > 0 it leads to the conclusion that all states (T, S ) in the support 
of 7r are reachable. Hence M is </>-irreducible. It is now enough to invoke the well-known 
fact that 0-irreducibility and aperiodicity imply ergodicity in total variation norm, see for 
example (Theorem 4, Roberts and Rosenthal, 2004). ■ 


Remark 4 Condition K(S,S') > 0 is clearly satisfied by FFBS algorithm, because it is 
equivalent to independent sampling from p(S\T , V, Y). This condition is also satisfied by the 
particle Gibbs algorithm described in the next section. 

In the case of CTBN we can use its dependence structure by introducing a Gibbs sampler 
over nodes of the graph (V,£). The same idea was exploited in (El-Hay et ah, 2008; Rao 
and Teh, 2013). By (8), the full conditional distribution of node w given rest of the graph 
has the density 

p(x w |x_ w ,y)(xi7(^(0)|x_ w (0))p(x u ,||x paM ) H p(x u \\x Mu) )L(Y\x ). ( 12 ) 

u£ch(w) 

The density p(X w \\X pa ^) corresponds to a piecewise homogeneous Markov process. The 
expression n«Gch(«j) P(X U ||-T pa ( u )) can be treated as a part of likelihood, similarly as L(Y\X). 
Note that the conditional initial distribution v(X w {f))\ A_„,(0)) may be replaced in formula 
(12) by the joint initial distribution z/(X(0)), because these two quanities are proportional 
as functions of X w (0). The step which leaves p(X w \X_ w ,Y) as invariant measure can be 
realized by the general algorithm described above. The Gibbs sampler for CTBN is sum¬ 
marized below. 

Algorithm 3 Gibbs sampler for CTBN. 

for io6 V (in a deterministic or random order) do 

Simulate ( T w , S w ) using single step of MCMC algorithm targeting p(X w \X_ w , Y), with 
X- w fixed. 

end for 


Note that if observations of different nodes are independent i.e. L(Y\X) = n^ev L W (Y\X W ) 
then the full conditional distributions defined by (12) reduce to 

p(X w \X_ w , Y) oc v(X(0))p(X w \\X Mw ))L w (Y\X w ) p(X u ||A pa ( u) ) . 

u€ch(w) 

Hence within the Gibbs sampler, the step for node w need not involve evaluation of full 
likelihood. An immediate corollary from Proposition 3 is that the Gibbs sampler for CTBN 
is also ergodic. Note that in the step of adding virtual jumps we can choose the intensity 
parameters A or 6 globally but, more generally, we can define different intensities for every 
node w, say X w or 6 W . 
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Corollary 5 Assume that for every node w G V we have \ w > max s „ sPa („) Q v (s pa ^; s w ) 
in the case of uniformization or 9 W > 0 for homogeneous virtual jumps. Then the Gibbs 
sampler for CTBN (with either the particle Gibbs or FFBS used in sampling of a new 
skeleton) is ergodic. 

7. Particle Gibbs 

In this section we will use notation which is standard in the literature on sequential Monte 
Carlo (SMC). Let SQ : k = (so, si,..., s^). Consider a sequence of unnormalized densities 
7fc(so:fc) on increasing product spaces S k+1 for k = 1,... n. The corresponding normalized 
probability densities are 


nk(so : k) 


7k(s 0:k) 

Zk 


where Zj. = f Jk( s O:k)dso-.k is the normalizing constant. A special case is the so-called 
state-space model where n k (s 0:k ) = p{so-.k\yo-.k ) and 7fc(s 0 :fc) = p(s 0:k , yo-.k)- In the sequel we 
consider the state-space model because the distribution of skeleton S given times of jumps 
(true and virtual; (T,V)) fits in this scheme c.f. (11). 

Particle MCMC methods introduced by Andrieu et al. (2010) provide a general frame¬ 
work for constructing an MCMC kernel targeting 7r(so:fc) with transition rule based on SMC 
algorithms. Before we describe particle MCMC in detail, we first have to recall standard 
SMC methods, e.g. Doucet et al. (2001); Doucet and Johansen (2009); Del Moral et al. 
(2006); Pitt and Shephard (1999). SMC sequentially approximates each of probability dis¬ 
tributions 7Tfc by an empirical distribution 

N i 

^k (dS0:fe) = J Vk 3 6 e 0 k ( dS0:fc ) ’ ( 13 ) 

i=i w i 

where {£o-fc, I s a weighted particle system. The system is propagated as follows. 

Given previous approximation {£o-fc-i’ w k-\}f=i a I b me k — 1, first we draw {£o fc-i}^=i 
from the multinomial distribution with probabilities proportional to weights {w^._ 1 }^ =1 
(this resampling step is introduced to avoid degeneracy of weights). Next we generate 
new particles {Cfc}£i according to r k (-l^o-.k-i) and set Q)-.k = (^o-.k-i^k)- Here r k is an 
instrumental kernel from S k to S (identified with a conditional density). Finally we compute 
new weights 


_7fc(£p : fc)_ 

7fe-l(^0:fe-l) r fe(^fello:fe-l) 


( 14 ) 


The SMC algorithm is summarized below. 
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Algorithm 4 SMC algorithm. 
Initialization 
for i = 1, ..., N do 

Draw Co ~ m(')- 
Compute weights 


end for 
Main loop 
for k = 1 ,... ,n do 
Resampling step: 

Draw ancestors 


7o(£q 

r o(C 0 ) 


Multinomial 



w 


k-l 


i 


and set 


£1 _ £ a t 

sO:fc-l — S0:fc-1 * 




Propagation step: 
for i = 1,..., N do 

Draw 

Ck ~ r k(-\C:k-l) 

and set 

Co-.k = (io:k-l’€l) ■ 


Compute weights w l k from (14). 

end for 
end for 


The particle Gibbs (PGS) algorithm introduced by Andrieu et al. (2010) is based on 
SMC algorithm. Given a previous path so:n we run an SMC algorithm with one path fixed, 
say Cy.n = s 0 :n, and obtain system of particles {^Q :n , w l n }^ =1 . Next we draw a new path s' 0 . n 
with probability 

P (4n = Co :n) « K ■ 

This procedure defines the following Markov kernel: 


K(S 0:n ,-)=E(^(-)\s k:n = ^ n ), 

where is defined by (13). It is shown that for every number of particles N larger than 
one, 7 r n is a stationary distribution for kernel K (Andrieu et al., 2010). 

There is a well-known problem of path degeneracy of SMC samplers (Doucet and Jo¬ 
hansen, 2009). For large n, the beginning of the path Sq.^ can be based on only few 
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trajectories. This may lead to poor mixing of the particle Gibbs sampler, because kernel 
K with high probability leaves the beginning of trajectory unchanged, see (Chopin and 
Singh, to appear 2015; Lindsten and Schon, 2013; Lindsten et ah, 2014b). A remedy for 
this problem can be an additional step of ancestor resampling proposed by Lindsten et al. 
(2014b). Let Cy n = so-. n be the fixed trajectory in the particle Gibbs algorithm. For ev¬ 
ery k = 1,... ,n we sample an ancestor of £ k from the set of trajectories {£o-fc-i}^=L with 
probabilities proportional to weights 


w 


k—l\n = w k-r 


7n((?0:fc-l’ s k:n)) 


(15) 


The above modification does not change the invariant measure of the Markov kernel (The¬ 
orem 1, Lindsten et al., 2014b). 

Now we are ready to describe precisely the step of sampling a new skeleton in the 
MCMC algorithm for hidden Markov jump processes. Let us recall (11). The conditional 
distribution of skeleton given moments of jumps and evidence is of the form 


p(S\T,V,Y) = v(s 0 )go(s 0 ) p ( s k-i, s k )g k (s k ) , 

k= 1 

where P is a transition matrix of some Markov chain. A standard choice is to use priors as 
instrumental kernels r i.e. tq = v and r k = P for k > 0. This choice leads to a simplified 
form of weights: 

w l o = 9o(Co) , K = g k {C k ) , w l k -\\ n = w k -iP{C k -i, s k ) , 

for i = 1..., N and k = 1,... ,n. If the algorithm is applied in a CTBN setting within 
a Gibbs sampler step, then the initial conditional distribution zz(A„,(0)|AT^O)) in (12) 
might be difficult to sample from. Then we can use a different instrumental distribution 
ro and compute weights w l 0 = <?o(£o) z/ (£o)/ r o(£o)- However, in many scenarios the initial 
configuration is deterministic (u is concentrated at a single configuration) and then this 
problem disappears. In algorithm 5 we summarize the particle Gibbs with ancestor sampling 
(PGAS). 

Remark 6 There exists also a Metropolis type of particle MCMC algorithm (PMH). How¬ 
ever, straightforward application of PMH within MCMC algorithm for Markov jump pro¬ 
cesses is not possible, because the acceptance probability in PMH depends on estimates of 
the normning constant Z n obtained in two consecutive steps. In the case of Markov jump 
processes, the dimension of the state space (n + 1) can be different in every step, because 
virtual jumps are either added or removed. Therefore to implement PMH in this context, 
transdimensional Metropolis type moves would be needed. 
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Algorithm 5 PGAS for sampling the new skeleton. 
Input: Current skeleton so :n . 

Output: New skeleton s' 0 . n . 

Set £():n = S0:n- 

Compute weights w k = gk(sk) for k = 0, 
for i = 1,..., N — 1 do 
Draw e 0 ~ *(■)■ 

Compute corresponding weights w l 0 = go(£o)- 

end for 

for k = 1,..., n do 
(Resampling) 

Draw ancestors 


and set 


Kh=i Multinomial I N, 


w 


k -1 


EM k -1 


<-i \ 


# _ fi a k 

S 0 :fc -1 “ S 0 :fe -1 ’ 


(Ancestor resampling) 

Draw J with probability 


and set 


P(J = i) 


w k-i p (?k-n s k) 


£():n — (£o:fc—1> s k:n) 


(Propagation step) 
for i = 1,..., N — 1 do 
Draw C k ~ ) 

Compute weights wf 

end for 
end for 

Draw I with probability 


and set = (Co : fc-n€fc)- 
= #*;(££)• 


P(7 = i) 



and set 
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8. Numerical experiments 

In this section we present results of simulations which demonstrate the efficiency of the 
proposed algorithm. We concentrate on the case of CTBNs. We start with a toy example 
with only two nodes joined by a single arrow (X —» Y), i.e. the simplest hidden Markov 
model with continuous time. Next we consider a CTBN with a chain structure. The 
last example is the Lotka-Volterra predator-prey model. To the best of our knowledge the 
algorithm of Rao and Teh (2013) is the most efficient method to deal with hidden continuous 
time Markov processes. For this reason we use it in our comparisons. In the simulations 
we use RCPP (Eddelbuettel et ah, 2011) implementation of algorithms. All our codes are 
available at a request. 


8.1 Toy example 

Consider a CTBN presented in Figure 1. Both nodes (X, Y) have two possible states, 
{1,2}. Node Y is fully observed, X is hidden (apart from the beginning and the end of its 
trajectory). The transition intensities of X are independent of Y and given by 


Qx 


-10 10 \ 
10 -10 ) ' 


The conditional intensities of Y are given by 


Qy\x=i 


-10 10 
10 -10 


Qy\X=2 


-100 100 \ 
100 -100 ) ■ 


Assume that we observe the full trajectory of Y in time interval [0,1] and we also observe 
X at moments 0 and 1. Our goal is to sample from the posterior distribution of X([0,1]) 
given the evidence (A(0), A'(l), Y([0,1])). We run our particle MCMC algorithm in two 
versions. Times of virtual jumps are added via uniformization with A = 20 in the first ver¬ 
sion and distributed according to the homogeneous Poisson process with intensity 9 = 10 
in the second version. In both cases the expected number of virtual jumps is approximately 
the same. We run our algorithm with PGAS with 2 and 4 particles. For a comparison we 
apply Rao and Teh (2013) algorithm with FFBS subroutine. The actual (hidden) trajec¬ 
tory A1([0,1]), the posterior means and standard deviation of MCMC approximation are 
presented in Figure 2. These results are based on 100 replications, each MCMC run of 
length 1000 with a burn-in time 100. We observe that the estimated trajectory is very sim¬ 
ilar for all the compared methods. The difference is in the variance and consequently in the 
width of the “confidence bands”. In this example the algorithm with uniformization outper¬ 
forms the algorithm with homogeneous virtual jumps. As it was expected, our algorithms 
with PGAS have the variance greater than those with FFBS but the difference between 
PGAS with 4 particles and FFBS is not too big. To analyze the rate of convergence of 
the algorithms we also compute standard deviation of sufficient statistics (number of jumps 
and time spent in each state). The results are shown in Figure 3. Again we obtain similar 
conclusions. It is clearly visible that the algorithms with FFBS have lower variance than 
those with PGAS. However, the difference is rapidly decreasing with increasing number of 
particles. Approximately the cost of FFBS is of order 0(\ ff| 2 E(n)) and the cost of PGAS 
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Figure 1: Toy example 


Uniformization 


B PGAS with 2 particles 


- 


PGAS with 4 particles - 


FFBS 



Virtual jumps from the homogenous Poisson process 

- PGAS with 2 particles- PGAS with 4 particles - FFBS 



Figure 2: Posterior mean and standard deviation of MCMC approximation for the toy ex¬ 
ample. 
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Standard deviation of estimates of sufficient statistics 



Number of particles 


- 2 . 


- FFBS 


Figure 3: Standard deviation of sufficient statistics versus number of iterations of MCMC 
algorithm. 


is 0(IVE(n)). Hence in the example under consideration, we obtain comparable quality of 
estimation for the algorithms with FFBS and PGAS which have the same computational 
cost. 

8.2 Chain network 

The model considered in this subsection is based on (Rao and Teh, 2013, Subsection 5.5). It 
is a network which consists of M nodes equipped with the “chain” graph 1 —^ 2 —> • • • —» M. 
For each node the set of possible states is {1,..., S'}. The transition intensities of every node, 
except the first one, depend on current state of the previous node. Namely, off-diagonal 
elements of intensity matrix are given by 

Q i(xi,x'i) 


Qm (%m— 1 j X m 

In words, the head node leaves the current state x\ with intensity 1 and prefers to jump to 
x\ T 1 mod S. For any other node, if its current state x m differs from that of parent state 
x m -\ then intensity of jumping is 2 and the process prefers to jump to x m -\. If x m = x m -\ 


= 


\ if x' x = (x\ + 1) mod S ; 

2 (g- 2 ) otherwise , 

] if x m — x m —i , 

1 if x m -i / x m and x' m = x m -i ; 
^32 otherwise. 
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Figure 4: Ratio of CPU time (FFBS/PGAS) needed to generate a sample with ESS = 100 
for the chain model: with increasing length of time interval (left), increasing 
number of nodes (center), increasing size of the state space (right). 


then the process leaves the current state with intensity 1 and chooses a new state at random. 
We observe the process at the beginning and at the end of time interval [0,T]. 

We run two MCMC algorithms targeting the posterior distribution over latent trajec¬ 
tories. Our algorithm with PGAS based on N = 10 particles is compared with Rao and 
Teh’s algorithm with FFBS. Both the algorithms use uniformization, with A equal to twice 
the intensity of leaving the current state. 

We begin with a network with M = 3 nodes, S = 2 states of every node and time 
length T = 5. In the first series of our experiments, we increase the number of nodes to 
M = 6,12, 24. In the second series we increase the length of time interval to T = 10, 20,40. 
Finally we change the size of the state space to S = 10, 50,100. For each of the scenarios we 
run 20 replications of each of the MCMC algorithms. We use CODA R package (Plummer 
et al., 2006) to estimate the effective sample size of the following statistics: time spent in 
each state and number of jumps, for nodes m = 1,..., M. The median of all these quantities 
serves as an estimator of ESS for the whole CTBN. In Figure 4 we present the ratio of CPU 
time needed to generate a sample with ESS equal to 100. In the beginning, when state 
space is of size S = 2, FFBS is more effective, around 4-5 times as fast as PGAS. This is 
what should be expected, since the cost of PGAS with 10 particles is higher than the cost 
of FFBS for a small space size. Also as expected, the ratio does not change significantly if 
the number of nodes increases. The same seems to happen if we increase the length of time 
interval. This last fact is slightly surprising, because the rate of convergence of particle 
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Figure 5: Predator-Prey model as CTBN. 


Gibbs depends on the number of jumps, see Andrieu et al. (2013); Lindsten et al. (2014a). 
However, if the size of the state space increases then the cost of our proposed method 
becomess significantly lower than that of Rao and Teh (2013) approach. For instance if 
S' = 50 then our algorithm is twice as fast as Rao and Teh (2013) and for S = 100 it is 
more than 9 times as fast. Experiments with a different number of particles (N = 5, 20, 50) 
lead to the same conclusions. 


8.3 Lotka-Volterra model 

The last example is Lotka-Volterra model (Wilkinson, 2009; Opper and Sanguinetti, 2008), 
which describes evolution of two interacting populations of prey and predator species. This 
process can be viewed as the two node CTBN shown in Figure 5. Let x and y represent the 
size of prey and predator populations, respectively. The transition intensities are given by 


Q({x, y}, {z + 1 , y}) = ax , Q({x, y}, {x - 1 , y}) = f3xy , 

Q({x, y}, {x, y + 1 }) = 5xy , Q({x, y}, {x, y - 1}) = jy , 


All other intensities are 0. The state space is infinite: {0,1,...} x {0,1,...}. Following Rao 
and Teh (2013) we set the parameters as follows: a = 7 = 5 x 10 ~ 4 and (3 = 5 = lx 10 -4 . 
We condsider the process in time interval [0, 3000] with known initial position and noisy 
observations Y(t) at discrete times uniformly spaced in interval [0,1500] with the likelihood 
given by 


p(Y(t)\X(t)) oc 


2 \X(t)-Y(t)\ + 1Q —6 
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Figure 6: Resuts of MCMC aproximation for Predator-Prey model. Black line - true path, 
blue line - posterior mean, shadow - 90% credible interval, red points - observa¬ 
tions. 


Given this evidence, we estimate the posterior distribution over sample paths of X using 
our sampler. We compute estimates of X in the interval [0,1500] (where observations Y are 
available) and also predict values of X in the interval [1500,3000]. Due to unboudedness 
of intensities we are not able to use uniformization technique and so we use homogenous 
virtual jumps with 8 = 30. For this choice of 6, the number of virtual jumps is approximately 
equal to the number of true jumps. We run MCMC simulation of length 1000 with 100 
initial iterations treated as burn-in time. In the PGAS we use 100 particles. The results of 
simulation are given in Figure 6. We conclude that the quality of estimates is the same as 
in Rao and Teh (2013). Note that in opposition to Rao and Teh’s algorithm we need not 
truncate the state space and the computational cost of our algorithm is significantly lower, 
namely C(100E(n)) for our method and (P(200 2 E(n)) for Rao and Teh, with state space 
truncted to {0,..., 200} x {0,..., 200} as suggested by these authors. 

9. Conclusions 

In the present paper we propose a new MCMC algorithm for sampling from the posterior 
distribution of hidden trajectory of a Markov jump process. The general idea is the same as 
in Rao and Teh (2013), namely we alternately add virtual jumps and update the skeleton of 
the process. The main novelty is that our algorithm uses PGAS to sample the skeleton, while 
Rao and Teh use FFBS. Thus instead of sampling exactly from a conditional distribution, 
we make a step of a Markov chain which preserves this distribution. This modification 
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has some disadvantages, as slower mixing of the entire procedure. However, there are also 
important advantages of our approach. Unlike previous methods our algorithm can be 
implemented even if the state space is infinite. The cost of a single step of the proposed 
algorithm does not depend on the size of the state space. Consequently, we can recommend 
our algorithm for problems where the space is either infinite or finite but large. If the size 
of the state space is not big, then our algorithm with PGAS converges slightly slower than 
the algorithm with FFBS. However, the difference between them is rapidly diminishing if 
we increase the number of particles in PGAS. 

In the present paper we describe an algorithm for time homogeneous processes, only to 
avoid too many technical details. However, the generalization to non-homogeneous processes 
is rather straightforward. For details of adding virtual jumps in non-homogeneous case 
we refer to Rao and Teh (2012). Since PGAS can easly deal with general state spaces, 
our algorithm can also be applied to piecewise deterministic Markov jump processes on 
general state spaces (i.e. processes which evolve deterministically between jumps and move 
according to some Markov kernel at moments of jumps). 

We note that an important issue is to find an optimal number of virtual jumps. Small 
number of virtual jumps can lead to poor mixing of moments of jumps. In the case of PGAS, 
large number of virtual jumps not only increases computational cost, as it is in the case of 
FFBS, but may have a negative impact on mixing of the whole algorithm. It is because the 
convergence rate of particle Gibbs depends on the length of simulated trajectory. 

Appendix A. 

Proof [Proof of Proposition 1] We are to check that if (T, S) has the distribution given by (3) 
then (: Tj,Sj ) is distributed according to (1). First we compute the distribution of waiting 
time for the next true jump. Without loss of generality we can assume that the previous 
jump occurred at time 0 and X (0) = s. To get the first true jump we generate subsequent 
moments of potential jumps t\,... ,ti... such that U — U-\ are i.i.d. from Exp(R(s)). The 
candidate is accepted with probability Q(s)/R(s). Hence 






Q(s) exp{— Q(s)u}du = 1 — exp{— Q(s)t} . 
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We have obtained an expression which is exactly the c.d.f. of waiting time for the next jump 
of process with intensity matrix Q. To conclude the proof, it is enough to note that 

P(si = s'|s 0 = s) = P(si = s'|s 0 = s, s 0 ± Sl) 

_ Q(s,s')/R(s) _ Q(s,s') 

Q(s)/R(s) Q(s) 


Proof [Proof of Corollary 2] By construction of the thinning procedure and by Proposition 1 
we have 


P(Vj\X(tj- 1 ) = a,tj-i,tj) 


p(tj-l,tj\s) 


' /t(s)KI + 1 


exp {-(tj - tj-i)R(s)} 


Q(s ) exp {-(tj - tj-i)Q(s) 


( R ( s ) - QO)) 1 ^ 1 exp {-(tj - tj-MRis) - Q(s))} . 
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